J15221 Evanthia TNBC - CyTOF marker expression
Ludmila Danilova & Joseph Tandurella
2023-05-01
Create figures for the clinical paper
All data is coming from tables for the paper
# labels for time points
timeLab = c("Baseline", "Post Run-in", "Week 8")
# set colors
recistCol = c("Non-responder" = "red","Responder" = "blue")
cbrCol = c("CR,PR,SD" = "blue","PD" = "red")
subtypeCol = c("HR+" = "red", "TNBC" = "blue")
timeCol = setNames(c("black", 'grey45', 'grey75'), timeLab)
# set root folder
rootDir = "./"
# Load data patient data
dta_pts = read.csv(paste0(rootDir,"input_tables/table_patient_info.csv"), row.names = 1)
dta_pts$ID = rownames(dta_pts)
#==========================
# file with TILs and PDL1 data
cytof_dat = read.csv(paste0(rootDir,"input_tables/table_cytof_marker_expression.csv"))
#markers = c("CD40", "CD103", "CCR5", "CD141")
markers = colnames(cytof_dat)[4:ncol(cytof_dat)]
#===========================
# add patient info
cytof_dat = cbind(cytof_dat, dta_pts[cytof_dat$patientID,])
# take only baseline and time 1 values
df_BvT1 <- cytof_dat %>% filter(timePoint %in% c('Baseline','Post Run-in'))
# take only time 1 and time 2 values
df_T1vT2 <- cytof_dat %>% filter(timePoint %in% c('Week 8','Post Run-in'))P-values
Pairwise Wilcoxon test was used for each time point comparing responders vs non-responders. And paired Wilcoxon test was used for baseline vs post run-in and post run-in vs week 8. We took the differences between all time points and then compare those differences between responders vs non-responders to see if there is significant difference in changes after treatment
# Calculate p-values from Wilcoxon test: paired for time points and not paired for response
res = c()
# split by clusters
cytof_dat_byClust = split(cytof_dat, f = factor(cytof_dat$cluster))
for(j in names(cytof_dat_byClust))
{
#
d = cytof_dat_byClust[[j]]
for (k in markers)
{
# create per patient matrix with three time points in columns
# some values would be NA
mat <- data.frame(matrix(nrow = length(unique(d$patientID)),
ncol = 3, dimnames = list(unique(d$patientID),
levels(factor(d$timePoint)))), check.names = F)
# fill in matrix
for(i in levels(factor(d$timePoint)))
{
d1 = d %>% filter(timePoint == i)
mat[d1$patientID, i] = d1[,k]
}
# add response
mat$RECIST = cytof_dat[rownames(mat),"RECIST"]
# create treatment effect
# substract baseline from post run-in and week 8
mat$delta_T1_B = mat[,2]-mat[,1]
mat$delta_T2_B = mat[,3]-mat[,1]
# and week 8 from post run-in
mat$delta_T2_T1 = mat[,3]-mat[,2]
# running paired test
# Baseline vs post run-in
p1 = wilcox.test(mat[,1], mat[,2], paired = T)$p.value
# post run-in vs week 8
p2 = wilcox.test(mat[,2], mat[,3], paired = T)$p.value
# test by response
# for baseline
p3 = wilcox.test(Baseline ~ RECIST, data = mat)$p.value
# for post run-in
p4 = wilcox.test(mat[,2] ~ mat[,"RECIST"], data = mat)$p.value
# for week 8
p5 = wilcox.test(mat[,3] ~ mat[,"RECIST"], data = mat)$p.value
# compare differences in treatment effect in responders vs non-responders
p6 = wilcox.test(delta_T1_B ~ RECIST, data = mat)$p.value
p7 = wilcox.test(delta_T2_B ~ RECIST, data = mat)$p.value
p8 = wilcox.test(delta_T2_T1 ~ RECIST, data = mat)$p.value
# combine all p-values and add means
res = rbind(res, c(marker = k, cluster = j, BvsT1 = round(p1,3),
T1vsT2 = round(p2,3), B_RvsNR = round(p3,3),
T1_RvsNR = round(p4,3), T2_RvsNR = round(p5,3),
d_T1_B_RvsNR = round(p6,3),
d_T2_B_RvsNR = round(p7,3),
d_T2_T1_RvsNR = round(p8,3),
mean_d_T1_B_R = round(mean(mat %>% filter(RECIST == "Responder") %>% select(delta_T1_B) %>% unlist(), na.rm = T),3),
mean_d_T1_B_NR = round(mean(mat %>% filter(RECIST == "Non-responder") %>% select(delta_T1_B) %>% unlist(), na.rm = T),3),
mean_d_T2_B_R = round(mean(mat %>% filter(RECIST == "Responder") %>% select(delta_T2_B) %>% unlist(), na.rm = T),3),
mean_d_T2_B_NR = round(mean(mat %>% filter(RECIST == "Non-responder") %>% select(delta_T2_B) %>% unlist(), na.rm = T),3),
mean_d_T2_T1_R = round(mean(mat %>% filter(RECIST == "Responder") %>% select(delta_T2_T1) %>% unlist(), na.rm = T),3),
mean_d_T2_T1_NR = round(mean(mat %>% filter(RECIST == "Non-responder") %>% select(delta_T2_T1) %>% unlist(), na.rm = T),3)))
}
}
res = res[order(res[,'marker']),]
datatable(res)# make boxplots for markers/clusters that have p-value < 0.05 in d_T1_B_RvsNR
sets = data.frame(res)
sets = sets %>% filter(d_T1_B_RvsNR < 0.05 & cluster%in% c('B','CD4_T','CD8_T','DNT','DPT','NK_like',"cDC1"))
for (i in 1:nrow(sets))
{
k = sets[i,'cluster']
j = sets[i,'marker']
d = cytof_dat %>% filter(cluster == k)
d = CreateAllFacet(d,"RECIST")
df = cbind(protein = d[,j],
d[,c("patientID", "timePoint", "RECIST","Subtype", "facet") ] )
p = ggpaired(df , x="timePoint",y="protein", id = "patientID",
line.size = 0.4, title = paste0(j,". ", k),
line.color = "RECIST",palette = recistCol, xlab="Time") +
facet_wrap("facet") +
scale_x_discrete(labels=timeLab)
print(p)
}
# additionally, plot CD103 in cDC1 cluster
k = "cDC1"
j = "CD103"
d = cytof_dat %>% filter(cluster == k)
d = CreateAllFacet(d,"RECIST")
df = cbind(protein = d[,j],
d[,c("patientID", "timePoint", "RECIST","Subtype", "facet") ] )
p = ggpaired(df , x="timePoint",y="protein", id = "patientID",
line.size = 0.4, title = paste0(j,". ", k),
line.color = "RECIST",palette = recistCol, xlab="Time") +
facet_wrap("facet") +
scale_x_discrete(labels=timeLab)
print(p)
Baseline vs Post Run-in
#ggpaired(df_BvT1 %>% filter(cluster %in% c("GMDSC")),
for (i in markers)
{
p = ggpaired(df_BvT1, x="timePoint",y=i, id = "patientID", line.color = "RECIST", line.size = 1,
palette = recistCol, title = paste0(i,". Baseline vs Post Run-in"), xlab="Time", ylab="MMI") +
scale_x_discrete(labels=timeLab)+
facet_wrap("cluster", scales = 'free') +
theme(axis.text.x=element_text(angle = 45, hjust = 1))
print(p)
}
Post Run-in vs Week 8
for (i in markers)
{
p = ggpaired(df_T1vT2, x="timePoint",y=i, id = "patientID", line.color = "RECIST", line.size = 1,
palette = recistCol, title = paste0(i,". Post Run-in vs Week 8"), xlab="Time", ylab="MMI") +
# scale_x_discrete(labels=timeLab)+
facet_wrap("cluster", scales = 'free') +
theme(axis.text.x=element_text(angle = 45, hjust = 1))
print(p)
}
sessionInfo()## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] scales_1.2.1 gridExtra_2.3 DT_0.27 data.table_1.14.8
## [5] ggrepel_0.9.3 ggpubr_0.6.0 ggplot2_3.4.2 dplyr_1.1.2
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 xfun_0.39 bslib_0.4.2 purrr_1.0.1
## [5] carData_3.0-5 colorspace_2.1-0 vctrs_0.6.2 generics_0.1.3
## [9] htmltools_0.5.5 yaml_2.3.7 utf8_1.2.3 rlang_1.1.0
## [13] jquerylib_0.1.4 pillar_1.9.0 glue_1.6.2 withr_2.5.0
## [17] lifecycle_1.0.3 munsell_0.5.0 ggsignif_0.6.4 gtable_0.3.3
## [21] htmlwidgets_1.6.2 evaluate_0.20 labeling_0.4.2 knitr_1.42
## [25] fastmap_1.1.1 crosstalk_1.2.0 fansi_1.0.4 highr_0.10
## [29] broom_1.0.4 Rcpp_1.0.10 backports_1.4.1 cachem_1.0.7
## [33] jsonlite_1.8.4 abind_1.4-5 farver_2.1.1 digest_0.6.31
## [37] rstatix_0.7.2 grid_4.2.2 cli_3.6.1 tools_4.2.2
## [41] magrittr_2.0.3 sass_0.4.5 tibble_3.2.1 tidyr_1.3.0
## [45] car_3.1-2 pkgconfig_2.0.3 ellipsis_0.3.2 rmarkdown_2.21
## [49] rstudioapi_0.14 R6_2.5.1 compiler_4.2.2